load('hdAll.mat');
hd=hd+hd';
K=5;n=212;Iter=10;
% center=[0, 10, 20, 30, 40, 49]+1;% 51张图 mesh编号中心位置
% center=[0, 4, 7, 3, 4, 5]+1;% 51张图 mesh编号中心位置
% all=randperm(n); center=all(1:K);
space=floor(n/K);center=1:space:n;
% center=floor(linspace(1,212,K));
C=[];

% % 第一次
% label=zeros(n,1);
% for i=1:length(center)
%     label(center(i))=i;
% end
% 
% %E步
% for i=1:n
%     if label(i)==0
%         if i<center(1)
%             label(i)=1;
%             continue;
%         end
%         if i>center(K)
%             label(i)=K;
%             continue;
%         end
%         in=center(1:end-1)<i&center(2:end)>i;
%         ind=find(in==1);
%         class2=[ind,ind+1];
%         center2=center(class2);
%         dist=hd(i,center2);
%         [~,index]=min(dist);
%         label(i)=class2(index);
%     end
% end
% 
% cost=0;
% for i=1:n
%     cost=cost+hd(i,center(label(i)));
% end
% C=[C cost];



for ii=1:10
    label=zeros(n,1);
    for i=1:length(center)
        label(center(i))=i;
    end
    
    %E步
    for i=1:n
        if label(i)==0
            if i<center(1)
                label(i)=1;
                continue;
            end
            if i>center(K)
                label(i)=K;
                continue;
            end
            in=center(1:end-1)<i&center(2:end)>i;
            ind=find(in==1);
            class2=[ind,ind+1];
            center2=center(class2);
            dist=hd(i,center2);
            [~,index]=min(dist);
            label(i)=class2(index);
        end
    end
    
    %混分的情况例如 1 1 1 1 2 1 1 2 2 2 1 2 2 2 2
    %M步骤
    for i=1:length(center)
        a=find(label==i);
        center(i)=a(ceil(length(a)/2));
    end
    %     center=sort(center);
    
    cost=0;
    for i=1:n
        cost=cost+hd(i,center(label(i)));
    end
    C=[C cost];
%     center
end
C
% center
% tabulate(label)

% T=[1:212;center(label)]';
% csvwrite(['keyframe',num2str(K),'.csv'],T)